Alternating dynamic state in intrinsic Josephson-junction stacks self-generated by 

internal resonance 
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Intrinsic Josephson-junction stacks realized in high-temperature superconductors provide a very 
attractive base for developing coherent sources of electromagnetic radiation in the terahertz fre- 
quency range. A promising way to synchronize phase oscillations in all the junctions is to excite 
an internal cavity resonance. We demonstrate that this resonance promotes the formation of an 
alternating coherent state, in which the system spontaneously splits into two subsystems with dif- 
ferent phase-oscillation patterns. There is a static phase shift between the oscillations in the two 
subsystems which changes from to 2tt in a narrow region near the stack center. The oscillating 
electric and magnetic fields are almost homogeneous in all the junctions. The formation of this 
state promotes efficient pumping of the energy into the cavity resonance leading to strong resonance 
features in the current- volt age dependence. 
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High-temperature superconductors, such as 
Bi2Sr2CaCu20s (BSCCO), are composed of two- 
dimensional superconducting CuCVlayers coupled via 
the Josephson effect. 1 The large packing density of 
these intrinsic junctions makes these compounds very 
attractive for developing coherent generators of elec- 
tromagnetic [em] radiation based on the ac Josephson 
effect. Moreover, a large value of the superconducting 
gap allows to bring the operation frequency of potential 
devices into the practically important terahertz range. 
To develop a powerful source, the major challenge is 
to synchronize the oscillations of the superconducting 
phases in a large number of junctions. A very promising 
route is to excite an internal cavity resonance in finite- 
size samples (mesas) j 2 -* 3 - which can entrain oscillations in 
a very large number of junctions. The frequency of this 
so-called in-phase Fiske mode is set by the lateral size 
of the mesa. The experimental demonstration of this 
mechanism 2 has brought the quest for superconducting 
terahertz sources to a new level. 

In general, a mechanism of pumping energy into the 
cavity mode is a nontrivial issue. Homogeneous phase 
oscillations at zero magnetic field do not couple to the 
Fiske modes. Such coupling can be facilitated by intro- 
ducing an external modulation of the Josephson critical 
current density. 3 In this case the amplitudes of the gen- 
erated standing wave and of the produced radiation are 
proportional to the strength of modulation. 

In this Letter we explore an interesting alternative 
possibility. Numerically solving the dynamics equation 
for the Josephson-junction stacks, we found that near 
the resonance an inhomogeneous synchronized state is 
formed. In this state, the system spontaneously splits 
into two subsystems with different phase-oscillation pat- 
terns, formally corresponding to fiuxon-antifluxon oscil- 
lations. Inspired by numerics, we also succeeded to build 
such solution analytically. The phase oscillations in two 
subsystems have a static phase shift which has a soliton- 
shape coordinate dependence, changing from at one 
side to 2tt at other side. This change occurs within the 
narrow region near the center of the stack and the width 



of this region shrinks when approaching to the resonance. 
In-spite of the difference in the phase oscillation patterns 
for the two subsystems, the oscillating electric and mag- 
netic fields are almost identical in all the junctions. The 
formation of this state strongly enhances coupling to the 
resonance mode and promotes efficient pumping of en- 
ergy into the cavity resonance. Such state was also found 
recently by Lin and Hu^ 

The dynamic equations for the Josephson-junction 
stacks have been derived in Refs.[5| and have been used in 
numerous simulation studies^, mostly to study dynamics 
of the Josephson vortex lattice formed by a magnetic field 
applied along the layers. We present these equations in 
the form of coupled time evolution equations for reduced 
electric and magnetic fields, e n and h n , phase differences, 
(f n , and the in-plane phase gradients, fc n , which are very 
convenient for numerical implementation, 

de n /dr = -v c e n - g{u) sin (p n + dh n /du, (la) 

d(p n /dr = e n , (lb) 

dk n /dr = -(k n + h n - ft n -i) /"ab, (lc) 

h n = f (d(p n /du - k n+1 + k n ) . (Id) 

The units and definitions of parameters are summa- 
rized in Table 1 and its caption. These reduced equa- 
tions depend on three parameters, v c = 47rt c / (e c uj p ), 
v ah = 47rcr ab /(e c u;p>y 2 ), and £ = A/s, where a c and cr ab 
are components of the quasiparticle conductivity. We 
neglected the in-plane displacement current which would 
give a term ~ d 2 k n /dr 2 , because relevant frequencies are 
much smaller than the in-plane plasma frequency. 

We simulated a stack containing N junctions (1 <n< 
A/"), having a width of LXj (0 < u < L), and assuming that 
the dynamic state is homogeneous in the third direction. 
We study the voltage range corresponding to the Joseph- 
son frequencies close to the lowest in-phase resonance fre- 
quency uo i = tt£/L. The function g(u) = l—2r(u—L/2)/L 
in Eq. (QJi) describes a linear modulation of the Joseph- 
son current density, which provides coupling to this mode 
for c-axis homogeneous oscillations. 3 Our purpose is to 
probe the qualitative structure of the resonance states 
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and, for simplicity, we assume nonradiative boundary through the stack, and metallic contacts at the top and 
conditions at the edges, k n = 0, dip n /du — ^I/2£ 2 , the bottom, ho — fcjv+i = 0. 
for u = 0, L, where / is the transport current flowing 
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Table 1. Units of physical variables. Here co p is plasma frequency, Aj is the Josephson length, s is the interlayer 
period, A is the in-plane London penetration depth, 7 is the anisotropy factor, and jj is the Josephson current density. 



Variable 


time, r 


coordinate, u 


phase gradient, k n 


electric field, e n 


magnetic field, h n 


current density, j 




Unit 


1/Up 


Aj 


l/Aj 


®olu p / (2ttcs) 


$ /(2^ 7 A 2 ) 


jj 





0.2 
0.18 
0.16 
0.14 
0.12 \- 



N=50, L=25, l=100,v c =0.01, 


v ab =0.2 | I 
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FIG. 1: Simulated and theoretical CVDs in the vicinity of 
the resonance voltage for different states with V = (e n ) =00. 
Small circles show the dependence obtained for r = 0.4 and 
c-axis homogeneous initial state, which remains homogeneous 
with increasing current. The theoretical CVD for this stated 
is shown by grey dashed line. Small filled diamonds show the 
CVD for r — 0.4 and inhomogeneous state, which was self- 
generated when small n- dependent perturbation was added to 
the phase in the beginning of run for each current value. The 
open squares represents the CVD for r = 0, when the inho- 
mogeneous state has been used as the initial state at starting 
current. The inhomogeneous state is not sensitive to modu- 
lation. Dotted black line shows the theoretical curve for the 
inhomogeneous state based on Eqs. J6} and ([T3]) . 



FIG. 2: Snapshots of the phase and electric field configura- 
tions in the 8 bottom junctions for the same parameters as 
in Fig. [I] and j = 0.18 marked by the arrow.— The phases are 
reduced to the interval [— 7r,7r], so that jumps formally cor- 
respond to centers of the Josephson vortices. The dynamic 
state corresponds to alternating nucleation, motion and an- 
nihilation of vortices in the even junctions and ant ivort ices 
in the odd junctions. Vortex velocities near the edges much 
larger than near the center. In between 3 rd and 4 th config- 
urations jumplike annihilation of fluxons and nucleation of 
antifluxons take place at the right side. The lower plots show 
that the electric field is homogeneous in all junctions and 
has space/time dependence corresponding to the fundamen- 
tal mode. The annihilation/nucleation events correspond to 
maxima of electric fields at the edges. 



Figure [T] shows the current-voltage dependences 
(CVDs) obtained for representative system parameters 
listed in the plot and for two values of the modulation 
parameter r = and 0.4. The dependences have been ob- 
tained with increasing current. We observe a strong res- 
onance enhancement of the current due to the excitation 
of the internal cavity resonance. As our simulations do 
not take into account thermal noise, the emerging state 
is sensitive to the initial configuration. If we start with 
a c-axis homogeneous state it remains homogeneous up 
to certain current. In this case the resonance is excited 
due to the finite modulation and it is well described by 
theory developed in Ref. 3. However, if we add a small n- 
dependent perturbation to the phase at the start of every 
current run, the homogeneous state blows up and the sys- 
tem organizes itself into a coherent inhomogeneous state. 
We also studied a system without modulation using an 
inhomogeneous state as initial state and found that the 



corresponding CVD is practically undistinguishable from 
the one for the modulated system. Therefore, the mod- 
ulation of the critical current density triggers the transi- 
tion to the inhomogeneous state but once being formed, 
this state is not sensitive to the modulation any more. 

To understand the nature of the inhomogeneous state, 
we show in Fig. [2] the time evolution of the phase and 
electric field for the 8 bottom junctions for j = 0.18. 
We see that the system splits into two alternating sub- 
systems with different phase dynamics, corresponding to 
fluxon/antifluxon oscillations. In the first half period vor- 
tices nucleate at the left side in even junctions, rapidly 
move to the center, then, after slow motion near the cen- 
ter, rapidly annihilate at the right side. Immediately af- 
ter that, in the second half period, antivortices nucleate 
at the right side in the odd junctions, move to the left 
in a similar way, and annihilate at the left side. In spite 
of the difference in the phase dynamics between the two 
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FIG. 3: Cosine of half static phase shift between phases in two 
subsystems at different currents. One can see that the phase 
shift has shape of soliton and it narrows with approaching 
to the resonance. Solid line shows theoretical curve based 
on Eq. (12]) and using Eqs. ©, (TO]), and (13]) for j = 0.13. 
This function determines the coupling of the homogeneous 
Josephson oscillations to the internal resonance mode. 



subsystems, the electric and magnetic fields are almost 
identical in all junctions. For the electric field, this can 
be seen from the lower plots of Fig. [2] The dominating 
contribution to the oscillating electric field is given by 
the fundamental cavity mode, e n ex cos(ttu/L). 

The homogeneity of the electric field implies that there 
is a static phase shift between the phase oscillations in 
the two subsystems. Figure [3] shows the cosine of half of 
this phase shift at different currents. One can see that 
the phase shift has the shape of a soliton and its width 
shrinks with increasing current. 

To study the dynamic state analytically, we assume 
that the system is split into two alternating subsystems, 
(/?2 TO +i = <£i, ^2™ = Introducing new variables, 

= ((pi + ty?2)/2 and tp- = if2 — (fi, and excluding 
other variables, we derive from Eqs. (1]) for £ ^> 1 



dr 2 
dr 2 



z dr 
dip- 
z ~d7 



du 2 
1 



-2 sin (y?_/2) cos ip+. 



• sin cp+ cos (cp-/2) , (2a) 
d \ d 2 (f- 

(2b) 



dr J du 2 



We now obtain a self-consistent approximate solution of 
these equations for the dynamic state when the Joseph- 
son frequency, uo = (e n ), is close to the resonance fre- 
quency, uj\ = tt£/L. We will show that cp- is almost 
static. In this case the equation for ip+ coincides with 
the phase equation for the Josephson junction with mod- 
ulated Josephson current density 3 with modulation func- 
tion g(u) = cos((/?_/2), see Fig. [3] Near the resonance 
frequency, we make the mode projection for 

ip+(u,r) = ojt + Re[^exp(— iuor)) cos(ttu / L) , (3) 

and, assuming |^| <C 1, we obtain 



with g_ = 



IV C UJ 



cos(ttu/L) cos (y?_/2) du 



(4) 
(5) 



being the coupling parameter. This solution determines 
the CVD, which takes into account resonance enhance- 
ment of the Josephson current 3 , 



j(V) « v c V + 



g 2 _v c V/i 



(a ,2 _ V 2y + {VcV) 



2 ' 



(6) 



To evaluate (p-(u, r), we split it into static and dy- 
namic parts, up-{u,r) = (p-(u) + (p-{u, r). Further anal- 
ysis shows that <p-(u,r) <C 1. The static part is deter- 
mined by 

d 2 up-jdu 2 - 8C+(u) sin(<p_/2) = (7) 
with C+(u) = (cos <p+). Using Eqs. (3]) and (4]), we obtain 

(8) 



C+(u) ~ C\ cos(7ru/L), 

withCl = _^M_ 9- 



2 (o,2 - W J) 2 + (l/cW) 2 



• (9) 



Consider the region near the midpoint, u = L/2, where 
the static cosine can be approximated by the linear func- 
tion C+(u) ~ Ci(tt/L)(L/2 — u). Using the substitution 

v = (u- L/2)/l 8 with l s = [L/ (8ttCi)] 1/3 , (10) 

we can reduce Eq. (7j) near the midpoint to the dimen- 
sionless form 



d 2 cp-/dv 2 + vsin(<p_/2) = 0. 



(11) 



This equation allows for soliton solution in which cp_ 
changes from to 2tt within \v\ ~ 1, corresponding to 
\u — L/2\ ~ l s . In the case l s <C L the linear expan- 
sion for C+(u) is valid in the soliton core and Eq. ([TT]) 
determines its shape with a very good accuracy. In the 
range uo 2 — uo 2 ^> v c uj the condition l s <C L is equiva- 
lent to uf(l -(uj/ujx) 2 ) <4tt 2 ^ 2 . As typically £ w 150, 
this condition is always satisfied in the interesting fre- 
quency range. The soliton solution has the symmetry 
(f-(v) — 2tt — cp-(-v). Numerically solving Eq. (TT]) , we 
can interpolate the solution as 



(p- (v) « 7r exp 



{-(V2/3) " 



\v\ + C v 



,3/2 



cl< 2 



} (12) 



for v<0 with C v « 0.5129. 

To evaluate the coupling constant (5]), we note that 
cos(y?_/2) changes from 1 to —1 within a narrow re- 
gion near the midpoint, — L/2\ ~l s meaning that, up 
to terms ~ (/ S /L) 2 , cos((/?_/2) can be approximated by 
sign(L/2— iz) which gives g_ ^4/7r. Surprisingly, this self- 
generated steplike modulation provides the maximum 
possible coupling to the resonance mode. Correction to 
g- due to the finite soliton width can be evaluated as 



4 f L ^ 2 /7TU\ / LP- \ 

5g- = -J cos {^) V cos 2 ~ l ) du ~ ~ - 464 



4^ 
L 2 



4 



Adding this correction and using Eqs. (fTQ|) and ([9]), the 
total coupling parameter can be written as 




1.817 



L 2 



-2/3 > 



(13) 



As g- ~ 1, the used linear approximation |^| <C 1 is valid 
only at \u 2 — co 2 \ > 1. 

To evaluate the time-dependent part of (/?_, we 
represent it in the complex form, (p_ (r, u) = 
Re [(p-(u) exp(— icjr)], and separating the time- 
dependent part of Eq. (|2b|h we derive the equation 
for the complex amplitude 



(w +iv c u;)ip_ + 



2 sin — 



In the range uo ^> 1 , we estimate (p- « 2 sin ((/?_ /2) /a; 2 <C 
1. The small amplitude of </?-(r, justifies usage of the 
static approximation for (/?_ in the equation for ip+. 

To compare analytical results with numerics, we show 
in Fig. [3] the theoretical result for cos(<^_/2) at j = 0.13 
based on Eq. ([12]) . It accurately describes the numerical 
data. The theoretical prediction for the CVD based on 
Eqs. (|6]) and ([T3]) is shown in Fig. [TJ The linear approxi- 
mation describes well the numerical data for voltages not 
too close to the resonance. Due to the enhancement of 
nonlinearities in the vicinity of the resonance, the ana- 
lytical result overestimates the current increase. 

The found state looks similar to the fluxon-antifluxon 
oscillations in a single junction^ These oscillations ap- 
pear as a result of a parametric instability of the homoge- 
neous oscillations 1 - and lead to so-called zero-field steps 
in the CVDs*ii The Josephson frequency of such step is 
twice the frequency of the involved resonance. In spite 
of the apparent similarity, there are essential qualitative 
differences. In the case of a stack, the frequency coin- 



cides with the resonance frequency. The dynamic config- 
urations are also very different. In the case of a single 
junction, a well-developed fluxon nucleates at one side, 
moves with Swihart velocity to the other side, converts to 
the antifluxon there, which then moves back again with 
constant velocity 9 . In our case, there is a region stati- 
cally located near the center where rapid phase change 
±7r is localized, corresponding to 2tt phase change of (/>_. 
As a consequence, the centers of fluxons and antifluxons, 
formally defined as points where the phases are commen- 
surate with ±7r, spend most time near the center and very 
rapidly jump to and from the edges, see Fig.[2j Moreover, 
the fluxon interpretation of our oscillations is somewhat 
artificial, as there are no well-defined localized soliton 
excitations moving across the junctions. 

The alternating state is a plausible candidate for the 
coherent state responsible for resonant terahertz emission 
observed in Ref. Q. Even though we did not use bound- 
ary conditions accounting for the radiation, it is clear 
that the generation of such a state would lead to pow- 
erful emission. In fact, the radiation does not influence 
much the structure of the internal states for short mesas, 
N < 100, and em emission can be approximately com- 
puted from the oscillating electric fields at the edges. 3 
For taller mesas, the radiation may contribute to the res- 
onance damping but we do not expect that it will destroy 
the coherent state. The experimental resonance features 
in the CVDs are much weaker then the theoretical ones. 
The possible mechanisms reducing the amplitude of the 
resonance include noise, c-axis inhomogeneities, and ad- 
ditional damping channels not taken into account by the 
theoretical model. 
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